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Abstract 

Integer-programming formulations for the design of 
symmetric and balanced laminated plates under biax- 
ial compression are presented. Both maximization of 
buckling load for given total thickness, and the mini- 
mization of total thickness subject to a buckling con- 
straint are formulated. The design variables that define 
the stacking sequence of the laminate are zero-one in- 
tegers. It is shown that the formulation results in a lin- 
ear optimization problem that can be solved on readily 
available software. This is in contrast to the continu- 
ous case, where the design variables are the thicknesses 
of layers with specified ply orientations, and the op- 
timization problem is nonlinear. Constraints on the 
stacking sequence such as a limit on the number of 
contiguous plies of the same orientation and limits on 
in-plane stiffnesses are easily accomodated. Examples 
are presented for graphite-epoxy plates under uniaxial 
and biaxial compression using a commercial software 
package based on the branch-and-bound algorithm. 

Introduction 

The design of laminated plates for maximum buckling 
load has drawn much attention in recent years (e.g., 
[l]-[7]). Typically the design variables are either the 
ply orientations of the layers or the thicknesses of lay- 
ers assumed to have a given ply orientation. However, 
in many practical applications the ply orientations that 
may be used are limited to 0-deg, 90-deg and ±45-deg, 
and the thicknesses of the layers are limited to inte- 
ger multiples of the lamina thickness. This means that 
the basic design problem is to determine the stacking 
sequence of the composite laminate — a problem which 
calls for integer programming techniques. 

Integer programming techniques are often quite costly, 
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and for this reason there have been several attempts 
to use ad-hoc techniques in applications to structural 
optimization (e.g., [8], [9]). However, the laminate de- 
sign problem (when classical lamination theory is used) 
is simple enough to permit the use of standard inte- 
ger programming techniques. Thus Mesquita and Ka- 
mat [10] and Olsen and Vanderplaats [11] have applied 
the popular branch and bound technique to the op- 
timization of composite laminates with thickness and 
ply orientation design variables subject to frequency 
or strength constraints. In Reference 10 the method 
was applied directly to the nonlinear problem, while in 
Reference 11 the nonlinear problem was solved as a se- 
quence of linearized problems. A similar approach was 
used by John and Ramakrishnan [12] for the design of 
trusses using a discrete set of sections. 

The objective of the present work is to show that the 
stacking sequence design of a laminated plate for buck- 
ling can be formulated as a linear problem by using 
ply-orientation-identity design variables. Thus, widely 
available software for the solution of linear integer pro- 
gramming problems can be used. Both the maximiza- 
tion of buckling load for specified total thickness and 
the dual problem of minimizing total thickness for spec- 
ified loading are studied. 

Analysis and Optimization Formulation 

A simply supported laminated plate under biaxial com- 
pression is shown in Figure 1. The loads per unit length 
in the x and y directions are XN x and XN y , respectively, 
with X being an amplitude parameter. The laminate is 
assumed to be symmetric and composed of Q-deg, 90- 
deg and ±45-deg plies. Each ply has a constant thick- 
ness t. For most of the examples in this paper, the lam- 
inate is also assumed to be balanced (i.e., the number 
of 45-deg plies is equal to the number of -45-deg plies). 
The laminate is composed of N p plies with a total thick- 
ness of h = N p t. However, because in some situations 
we will not know the number of plies (this will be de- 
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termined by the optimization process) the number of 
plies is assumed to be smaller than an upper limit N. 
The laminate buckles when the load amplitude reaches 
a critical value A cr given as 


A (:r (r 7 j 1 n)/ 7 r 2 = 

Dll (?)* + mi + 2D66 ) (^) 2 (f ) 2 + ^22 (?)* 

(?) a *. + (*)’** 

(!) 

where m and n are the number of half waves in the x 
and y directions, respectively, selected so as to mini- 
mize A cr . In the present study the minimization over 
m and n is performed by checking for all value of m 
between 1 and m/, and all values of n between 1 and 
nj. The flexural stiffnesses Du, Dn, D 22 and D^e can 
be expressed in terms of three integrals, V 0 , Vi, and 
V 3 , and five material invariants t/,-, i = I, ...» 5, which 
depend on the stacking sequence [13] as 

Du = U\Vq + U 2 V\ + U 3 V 3 
D 22 — f/i Vq - D 2 Vi + U 3 V 3 
D l2 C/ 4 V 0 - D 3 V 3 

D 66 = D 5 V 0 - f/ 3 V 3 
where Vq, Vi and V3 are given as 
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It is convenient to work in terms of non dimensional 
loads n r , n y , flexural stiffnesses dij , integrals Do, tq 
and U 3 , and material constants u, defined as 
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Then A cr is given as 
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The nondirnensional flexural stiffnesses are given as 


d u — tljUo + t* 2 t>l + U 3 U 3 
d 2 2 - ** 1*>0 - ** 2*>1 + ** 3**3 
di 2 = U4V0 ~ U3V3 
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where h is the total thickness of the laminate, z is the 
distance from the plane of symmetry (see Figure 1), 6 
is the ply orientation angle, and pjt is a variable which 
is equal to one if the ith ply is occupied and is equal to 
zero if the ply is empty. Constraints are applied during 
the optimization to ensure that p* can be zero only for 
the outermost plies. The material invariants are 

U\ = -(3Qn -f 3Q 2 2 + 2Qi2 + 4<2 66 ) 

o 

f/ 2 = -{Qu - Q22) 

U3 — d" Qi 2 “ 2 Qi 2 — 4 Qee) (A) 

U4 — r(Qu + Q22 + 6Q12 " 4 Q$e) 

o 

f/5 = — (Qn + Q22 — 2Q12 + 4 Qse) 
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Because the laminate is symmetric only the plies be- 
low the plane of symmetry need to be defined. The 
ply stacking sequence is defined in terms of four sets 
of ply-orientation-identity variables o,*, n t , f? and f - n , 
i = 1, ■ • • , N/2 that are zero-one integer variables. The 
variable o,, n t , ff or /, m is equal to one if there is a 
0-deg, 90-deg, 45-deg or -45-deg ply, respectively, in 
the zth layer. Unlike conventional practice, it is more 
convenient here to number the plies so that the first 
one (i = 1 ) is nearest the plane of symmetry of the 
laminate, and the last one is on the outside (i = N/ 2 ). 
The stacking-sequence variables are used to express the 
nondirnensional integrals vq, and V 3 as 

Nl 2 

». = Ep»[(T) s - (2 r) 3 ] 
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where /£ and /™ do not appear in the expression for 
ri since the cosine of 90 degrees is equal to zero. Two 
optimization problems are formulated. The first is the 
optimization of a laminate with a fixed thickness for 
maximum buckling load, and the second is the opti- 
mization of a laminate with a given buckling load for 
minimum thickness. For the first optimization problem 
the lowest (over values of m and n) buckling load A* is 
maximized. The objective A* is not a smooth function 
of the design variables, and the standard device for re- 
moving this problem is to add A* as a design variable 
and require it to be less than or equal to each A cr (m, n). 
Thus, the optimization problem is formulated as 


find A*, and o it n;, ff //" 

i= 1, — ,A72 

to maximize A* 

such that A* < A cr (m,n), 

m- 1, n - 1 , ■ ■ • , n/ (12) 

<H + m + f[ + fr = 1 

i = l,",N/2 


N/2 

and - IT - 0 

1=1 

where the last constraint ensures that the number of 
45-deg and -45-deg plies is the same, so that the lami- 
nate is balanced. Equations (9)- (11) are used to calcu- 
late the nondimensiona! buckling load which is clearly 
a linear function of the stacking-sequence design vari- 
able. Therefore, the optimization problem (12) is a 
linear integer programming problem. 


This linear formulation should be contrasted to the one 
obtained when ply thicknesses are used as design vari- 
ables. That formulation results in a nonlinear optimiza- 
tion problem. Thus the continuous formulation results 
in a problem which is more difficult than the integer 
formulation presented here. 


specified buckling load, the number of plies is not spec- 
ified. However, the dual formulation is only slightly 
more complex. The number of plies, N } is selected 
large enough to insure that a laminate that does not 
buckle can be found. This can be done by analyzing 
a trial design and then scaling the laminate thickness 
so it does not buckle (the buckling load is proportional 
to the cube of the laminate thickness if the same ply 
stacking-sequence is repeated again and again). Now 
the laminate is designed permitting some of the outer 
layers to be empty. Buckling will not occur for a spec- 
ified N x and N v if A cr is greater than or equal to one. 
The problem is formulated as 


find o<, n t , ff rr 

i= 

N/2 

to minimize + n i + ff + fi n ) 

i=i 

such that A cr (m,n) > 1, 

i7i=l,- n 

Oi+ni + Jf + /r < 1 (13) 

i= l > .-.,JV/2 

N/2 

£(/?-r) = o 

»= i 

and Oi + + ff + ff 1 

< Oi - i + Hi-1 + ff- 1 + ffli 

i = 2 t .~,N/2 


where the last constraint ensures that if there are empty 
plies they are on the outside. 

In general, the solution to the optimization problem 
(13) is not unique. For example, the non-integer solu- 
tion could require 8.1 plies. The design from (13) will 
have 10 plies (N must be even because of symmetry), 
and it will have a substantial margin, that is A cr will be 
significantly larger than 1. Any weaker 10-ply design, 
that is one that has a A cr closer to 1, is also a legitimate 
solution of (13) in that it satisfies all the constraints 
and has the same value of the objective function. In 
the present work, to achieve a unique solution, it is as- 
sumed that the best design is the minimum thickness 
plate that has also the largest possible buckling mar- 
gin of all plates of the same thickness. To achieve this 
goal the objective function of (13) was modified by sub- 
stracting e A cr from it, where c is a small number (0.001 
for the results presented in the next section). 


Unlike the formulation of Eq. (12), in the dual prob- 
lem of minimizing the laminate thickness subject to a 
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Another reason for a nonunique solution is that in terms 
of the calculation of the flexural stiffnesses of Eq. (2) 



there is no difference between the contribution of 45-deg 
and -45-deg plies. However, the buckling-load calcula- 
tion, Eq. (1), is not accurate for large values of D\e and 
Z? 26 * These can be minimized by selecting the p ositio ns 
of the 45-deg and -45-deg plies so as to minimize their 
combined contribution to Vo, Eq. (3). This selection 
was done by modifying manually the optimum design. 

In some cases it may be desirable to impose constraints 
on the stiffness of the plate. In the present study a 
limit on the in-plane stiffness in the x direction A\\ 
was considered as an example of such constraints. A 
constraint requiring A n to have a minimum value of 
A° lx can be written as 

>4iiMn-l>Q ( 14 ) 

As shown in the Appendix, this constraint can be ex- 
pressed as a linear function of the ply identity design 
variables in a manner similar to the buckling constraint 
(in Eqs. (9)-(ll)). 

Results 

Results were obtained for graphite-epoxy laminates 
(£1 = 18.5 x 10 6 psi, E 2 = 1.89 x 10 6 psi, G u = 
0.93 x 10 6 psi, i/i 2 = 0.3, t - 0.005 in.). The compu- 
tations were performed with the LIN DO program [14] 
which employs the branch-and-bound algorithm. First 
uniaxial loading was applied, and the buckling load was 
maximized for various plate aspect ratios ( a/b ) for lam- 
inates with 16 plies. It is known (e.g. [15]) that for low 
aspect ratios the optimum ply angle is 0 -deg, while for 
a/b larger than about 0.7 the optimum ply orientation 
is close to ±45-deg. This can also be expected from Eq. 
(9) since for a/b larger than 0.7, dee is the most impor- 
tant stiffness coefficient. A check was performed to see 
whether there was a transition region where the opti- 
mum stacking sequence would include both 0 -deg and 
±45-deg plies. It was found that if such a transition re- 
gion exists it is extremely narrow, in that even changes 
in the fourth significant digit of the aspect ratio were 
not fine enough to locate it. When the number of plies 
(N) was not divisible by four, so that a balanced ±45- 
deg laminate was not possible, the optimizer placed two 
0 -deg plies near the plane of symmetry of the laminate, 
as expected (because these less efficient plies have the 
smallest effect on dee there). 

Next the biaxial loading case was solved, and the re- 
sults are presented in Figure 2. It is known (e.g., [15]) 
that for aspect ratios less than 1.5 the optimum ply 
orientation is the same as for the uniaxial case, and for 
aspect ratios greater than 1.5, the value of the opti- 
mum ply angle increases rapidly as N y /N x increases, 


and that for large N y /N x , the optimum ply angle is 
90 degrees. Therefore, the case of biaxial loading for a 
laminated plate with an aspect ratio of 2 was selected. 
The reference axial load N x was fixed at 1 lb/in. (so 
that A cr is equal to the critical value of N x ), and the 
reference transverse load was increased from 0.1 to 3.0 
lb/in. The plate was specified to have 16 plies. For 
this choice, the value of A cr is equal to the value of N x 
at buckling. Two transition regions were found: one 
for N y /N x between 0.125 and 0.15 and the other for 
N y /N x between 2.4 and 2.45. The first region marked 
the transition from all ±45-deg plies to a combination 
of 90-deg and ±45-deg plies. The second region marked 
the transition from a combination of 90-deg and ±45- 
deg plies to all 90-deg plies. As the ratio N y /N x in- 
creased to 0.15 first two 90-deg plies were added, then 
four 90-deg plies (Af y /N x =0.25), then six 90-deg plies 
( Ny j N x — 1 ) , and finally all 90-deg plies (A r y /A r x = 2.45). 
As the transverse load N y became larger than the axial 
load N Xi the ± 45 -deg plies moved closer to the plane 
of symmetry until only 90-deg plies were present. This 
behavior is expected because when N y dominates, the 
plate behaves like a plate of aspect ratio of 0.5 under 
uniaxial load, and for that case the optimum angle is 
in the direction of the loading. 

When the number of contiguous plies in the same direc- 
tion is large, composite laminates are known to experi- 
ence matrix cracking. Therefore, it is desirable to limit 
the number of such contiguous plies. To demonstrate 
that such constraints can be easily added to the present 
formulation we imposed this constraint on the design 
obtained for N y /N x - 2 which had 5 contiguous 90-deg 
plies. This was implemented by adding the constraint 

^4 + n 5 ± ^6 ± 717 ± ns < 4 (14) 

The designs with and without this constraint are com- 
pared in Figure 3. It is seen that the penalty for limiting 
the number of contiguous plies is quite small. 

The case of N y /N x — 2 was used also for the pur- 
pose of checking on other aspects of the optimization. 
The first was the effect of the balanced laminate re- 
quirement. When this requirement was removed, the 
optimization selected a design with three 45-deg and 
five 90-deg plies. However, the buckling load changed 
by less than one hundredth of one percent. A second 
aspect was the effect of requiring that the design vari- 
ables be integers. Noninteger design variables describe 
hybrid plies. For example, o x = 0.5, n x = 0.5 means 
that the first ply has properties which are the aver- 
age of the elastic properties of zero and ninety degree 
material. When the requirement that the ply-identity 
variables be integers was removed the solution included 
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two hybrid plies. For i = 1 the ply was 70 percent 45- 
deg and 30 percent 90-deg, and for i = 4 the ply was 
70 percent -45-deg and 30 percent 90-deg, with the re- 
maining plies being 90-deg. The effect on the buckling 
load was again quite small — less than five hundredths 
of one percent. 

Another aspect of the optimization checked for this case 
was the effect of introducing a minimum stiffness re- 
quirement. The optimum laminate for this case, being 
dominated by 90-deg plies, has only 16 percent of the 
axial stiffness An of an all 0-deg laminate. A require- 
ment that An is at least 50 percent of the unidirectional 
laminate was added, with and without the requirement 
of no more than four contiguous plies. The results are 
compared to the original design in Figure 4. It is seen 
that the stiffness requirement is satifsfied by putting 0- 
deg plies near the plane of symmetry where they have 
only a minimal effect on the bending stiffnesses, and 
hence on the buckling load. The reduction in the buck- 
ling load is about 8 percent. This time effect of adding 
the requirements of no more than 4 contiguous plies had 
a nontrivial effect (7 percent reduction) on the buckling 
load. 

Next we solved the minimum thickness problem for a 
laminate with the same dimension. The axial load, N Xi 
was fixed at 30 lb /in., and the transverse load N y was 
varied from zero to 75 lb/in. The results are summa- 
rized in Figure 5. For N y = 0 we have a 10-ply laminate 
which is dominated by ±45-deg plies, with two 0-deg 
plies near the plane of symmetry. As N y is increased, 
the number of plies increases, and the laminate becomes 
dominated by 90-deg plies. However, the requirement 
of a balanced laminate tends to disturb the progression 
toward increasing number of 90-deg plies. For example, 
with loads that result in 12-ply laminates we can have 
either 4 or 8 ±45-deg plies, and the optimizer chooses 4, 
because 8 would leave only four 90-deg plies. However, 
when we increase the load so that we require 14-ply 
laminates, the number of ±45-deg plies jumps from 4 
to 8, because we can have now six 90-deg plies. 

Concluding Remarks 

The problem of stacking sequence design of composite 
laminates for minimum thickness subject to a buck- 
ling constraint or maximum buckling load for a given 
thickness was addressed. It w T as shown that the use of 
ply-orientation-identity design variables results in a lin- 
ear formulation of the problem unlike the use of more 
traditional ply-thickness design variables which lead to 
nonlinear formulation. The linear integer-programming 
formulation was solved using a commercially available 


program based on the branch-and-bound algorithm. It 
was also shown that the formulation can accomodate 
constraints on stiffnessess as well as constraints on the 
maximum number of contiguous plies of same angle. 
Results were presented for both uniaxial and biaxial 
loadings. 
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and the nondimensional integrals can be expressed in 
terms of the ply-identity design variables as 

N/2 
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»M = 2£K-n k ) (Al) 

k=l 
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«3 A = 2 ^(°Jk + n k ~ /jfc _ /**) 

*= 1 

In the example used in the Results section the lower 
limit on An is a specified fraction / of the correspond- 
ing stiffness of an n-ply all 0° laminate. For such a 
laminate vqa = = v$ A = Af , so that the constraint 

of Eq. (14) becomes 


All 

N(Ui + «2 + « 3 ) 


-/>o 
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Appendix — In-Plane Stiffness Constraint 


This appendix shows how a limit on the in-plane stiff- 
ness Ah can be formulated as a linear function of 
the ply-orientation-identity design variables. Limits on 
other stiffness components can be formulated in a sim- 
ilar way. 


The in-plane stiffness An is given as 


where 


and 


Ail = U\Vqa + V 2 V\A + U 3 V 3 A 

(A 1) 
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/ dz = 2 Pk 
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Ny = 0.25 Ny=0.5 

(45% -45” 2 , 90" 2 , 45% 90° 2 ) s (45% 90% -45% 90 - 3 . -45% 45°), 

X cr = 127.49 X er = 94.29 



N y = 1.0 

(90% 45% 90 Y *45% 90° ), 
Xcr = 61.99 



Ny =1.5 

(90” 3 , 45% 90* 2 , -45% 90°), 
Xcr s 46.18 


Laminate 

1 plane of 
symmetry 

Z 

Figure 1. Laminated plate geometry and 
loading. 


NyS 0.125 N y = 0.15 

(45% *45*2 » 45° , -45°, 45* 2 , •45") # (45°, -45° 2 , 45* , -45’, 90°, 45% 90*), 

X cr = 154.06 X cr = 148.46 



Figure 2. Maximum buckling load results for 
graphite-epoxy plate (16 plies, 
a=20 In, b=10 In, N =1 Ib/ln). 






Ny = 2.0 

(90‘ 5 .45%-45%90») s 
X cr = 36.84 



Ny r 2.1 

(90° 5 , 45°, 90°, -45°) s 
X w = 35.40 



Figure 2. Concluded. 
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(90° 5 , 45°, -45°, 90°) s 
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Ny = 0.0 lb/ln 
(45°, -45° 2 , 45°, 0°) s 
X cr = 1-41 
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N y = 7.5 lb/ln 
(-45°, 45°, 90° 3 ), 
X cr = 1.03 
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N y = 15 lb/ln 
(•45°, 90° 2 , 45°, 90°, 0°), 
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Ny = 22.5 lb/ln 
(90°, -45°, 90°, 45°, 90° 2 ), 
Xcr “ 1-05 


Figure 3. Effect of constraint of no more 
than four contiguous plies in 
same direction on design for 
N y /N x =2. 
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Ny = 45 lb/ln 

(90*3,45°, 90°, -45°, 90°), 
X cf = 1.02 
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Ny = 60 lb/ln 
(90° 5> 45°, -45°, 90°), 
X cr = 1.22 


N y = 75 lb/ln 

(9°‘l 6 ) 

X cr = 1-02 
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Figure 4. Effect of stiffness requirement 
on laminate design for N y /N x =2, 

N y =2.0. 


Figure 5. Minimum thickness results for 
graphite-epoxy plate (a=20 In, 
b=10 In, N =30 lb/ln). 
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